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, We recently derived explicit solutions of the leading-order Dokshitzer-Gribov-Lipatov-Altarelli- 

^ ■ Parisi (DGLAP) equations for the evolution of the singlet structure function Fs{x,Q^) and the 

, O ' gluon distribution G{x, Q^) using very efficient Laplace transform techniques. We apply our results 

, here to a study of the HERA data on deep inelastic ep scattering as recently combined by the HI 

and ZEUS groups. We use initial distributions F2^{x,Qq) and G{x,Qq) determined for x < 0.1 
by a global fit to the HERA data, and extended to a; = 1 using the shapes of those distributions 
determined in the CTEQ6L and MSTW2008LO analyses from fits to other data. Our final results 
, are insensitive at small x to the details of the extension. We obtain the singlet quark distribution 

Qh' Fs{x, Qo) from F^^ix, Qg) using small nonsinglet quark distributions taken from either the CTEQ6L 

or the MSTW2008LO analyses, evolve Fs and G to arbitrary Q^, and then convert the results to 
individual quark distributions. Finally, we show directly from a study of systematic trends in a 
comparison of the evolved F2^(x,Q'^) with the HERA data, that the assumption of leading-order 
DGLAP evolution is inconsistent with those data. 
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^ , PACS numbers: 12.38.Bx,12.38.-t,13.60.Hb 
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^ . I. INTRODUCTION 

' In recent papers [ll, we showed that it is possible to solve the coupled leading-order (LO) Dokshitzer-Gribov- 
Lipatov-Altarelli-Parisi (DGLAP) evolution equations 0-0] for the singlet quark structure function Fs{x,Q'^) = 
J2i x[qi{x, Q^) + qi{x, Q^)] and the gluon distribution G{x, Q^) = xg{x, Q^) directly using a method based on Laplace 
transforms. While the method is formally equivalent through the known connection between Laplace and Mellin 
' transforms @ to methods based on the latter - see, e.g. [1, Q , we find the present approach to be clearer intuitively 
' and much more efficient numerically. In particular, the distributions Fs{x, Q^) and G{x, Q^) at a virtuality can be 
^ , expressed simply as convolutions of the distributions Fs{x, Qq) and G{x, Qo) at a starting value Qq with analytically 
defined kernels in the ordinary variables. Alternatively, the results can be expressed as inverse Laplace transforms of 
products of the kernels in Laplace space with the Laplace transforms of the initial distributions. 

We perform the inverse Laplace transforms necessary in our approach using very fast and accurate new numerical 
algorithms Q . These do not require that we work on a preassigned numerical grid, and make the solution of the 
evolution equations at arbitrary values x and straightforward on desktop or laptop computers. We have extended 
the Laplace method elsewhere [3 to next-to- leading order (NLO) in ag, including to nonsinglet distributions, but will 
not pursue that extension here. 

In the present paper, we apply these methods to test the consistency of the assumed LO evolution of the structure 
functions with the HERA data jl0l - [l2l | on deep inelastic ep (or 7* p) s cattering, using those data as recently combined 
by the HI and ZEUS experimental groups [l^. As shown earlier [l4l [TH. if a LO treatment of the DGLAP evolution 
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is sufficient, the necessary starting distribution Go(x) = G(x,Qq) can be obtained from a global fit to the structure 
function Q^) by requiring that the LO evolution equation for F^^{x,Q'^) be satisfied for = Qq. Both 

F2^{x, Qq) and G{x, Qq) are then determined directly by experiment. 

To obtain our starting distributions, we perform the required global fit to using the HERA data for x < 0.1, 
and extend the fit to a: = 1 using the shape of that distribution as determined in the CTEQ6L [1^ and MSTW2008LO 
(l7l | analyses which included other DIS data at large x. Our final results at small x are insensitive to the details of the 
extension. We pick as a starting value for the Q^ evolution a value Qq = 4.5 GeV^, which is well within the region of 
dense data, and determine the starting G as described above. 

The singlet distribution Fs{x, Q^) differs from F^^(a;, Q^) by small nonsinglet contributions that depend primarily 
on the valence quark distributions, which agree fairly well for different LO analyses at moderate Q^ (c omp are, e.g. 
[Hi and [l3|). We wiU therefore simply use the results of the CTEQ6L and MSTW2008LO analyses [lllH to make 
the necessary conversion from to Fg at Qq = 4.5 GeV^, and the evolved nonsinglet contributions to convert the 
evolved Fs{x, Q^) back to the function F2^{x, Q^) which can be compared to the HERA data for Q^ ^ Qq. 

We also combine the evolved Fg with the nonsinglet distributions of CTEQ6L and MSTW2008LO to obtain a new 
set of CTEQ6L-like or MSTW-like quark distributions. Even though we use the same nonsinglet distributions as those 
authors, our final results differ from the originals because of our use of the combined HERA data rather than the 
original HI and ZEUS results, and, importantly, our use of starting distributions Fs{x,Qq) and G(a;, Qq) determined 
directly from experiment up to the small nonsinglet contributions to the former. 

We find that the evolved Q^) calculated using LO DGLAP evolution differs systematically in its dependence 

on X and Q^ from the combined HERA data at values of Q^ away from Qq. We conclude that LO DGLAP evolution 
is not consistent with the data, a conclusion reached less directly by other authors, e.g., in emphasize 
in this connection that the only fitting involved in our approach is in the QCD-independent global fit to the data on 
i^^^; we do not need to solve the complete set of evolution equations and then attempt to fit the data using the many 
input parameters typically introduced in the parameterization of initial parton distributions. 

Our conclusion on the inconsistency of LO evolution is not surprising. Next-to-leading-order (NLO) effects on the 
evolution are known to be large. However, our results give a direct demonstration of the necessity of going beyond 
LO independent of the substantial complications that a NLO analysis entails. 

In the Appendix, we present an accurate alternative method of testing LO evolution based on the exact LO evolution 
equation for Q^), and an approximate evolution equation for G{x, Q^). Its advantage is that the input necessary 

to test the assumption of LO evolution can be obtained directly from the measured F2^{x,Q'^). The application of 
this method to the HERA data leads to the same conclusion as stated above: the assumption of LO evolution is 
inconsistent with HERA data. 



In the present paper, we use the method developed in detail in |2| to solve the coupled DGLAP evolution 
equations for Fg and G. We will not give the details here, but note that our method is based on Laplace transforms. 
We first rewrite the evolution equations in terms of the variables v = ln(l/a;) and Q^ instead of x and Q^. The 
integral coupling terms in the equations then reduce to a form that involves convolutions in v, and the equations 
can be converted by Laplace transformation to factored homogeneous first-order differential equations in Q^ and a 
Laplace variable s, and solved directly. 

Using the notation Fs{v, Q^) = Fs{e^^ , Q^), G(u, Q^) = G{e~'" , Q^) for the distributions written in terms of v and 
Q^, and introducing the Laplace transforms 



we find that the Laplace-space distributions generated by evolution from Qq to Q^ can be expressed in terms of the 
initial distributions /so(s) = /s(s, Qq) and go(s) = .g(s, Qq) as 



II. PRELIMINARIES 



A. Solution of the coupled evolution equations for Fs and G 



fs{s,Q^) = C Fs{v,Q^);s 



(1) 



The kernels fc(s,r) in Eqs, 
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variable 

r{Q\Ql)^^j asiQ'')d(lnQ'^), (4) 

which vanishes for = Qq, with kff{s,0) = kgg{s,0) = 1 and kfg{s,0) = fcg/(s,0) = 0. The kernels also depend on 
the number nf oi active quarks. 

If we have parametrized the initial distributions accurately analytically, and Laplace transformed the results to 
obtain /so(s) and go{s), we can calculate the inverse Laplace transforms of fsis,Q^) and g{s,Q^) directly to obtain 
the evolved distributions Fs{v,Q^) and G{v,Q^), with 

Fs{v,Q^) = C-'{[kff{s,T)U{s);v] + [kfg{s,T)go{s);v]}, (5) 
G(«,g2) = C~'{[kgf{s,T)U{s);v] + [kggis,T)go{s);v]}. (6) 

Alternatively, using the convolution theorem to write the transforms of the products on the right-hand sides as 
convolutions, and using the fact that the inverse transforms of fso{s) and go{s) are the initial v-space distributions 
Fsoiv) = Fs{v, Qq), Go{v) = G{v, Qq), we can write the solutions in the more intuitive form 



F,{v,Q^) = / KFF{v~w,T{Q\Ql))F,o{w)dw+ Kpc {v - w,t{Q\QI)) Go{w) dw, (7) 
Jo Jo 

G{v,Q^) ^ / KGF{v-w,T{Q^,Ql))F,o{w)dw+ Kgg {v - w,t{Q\qI)) G„{w) dw, (8) 



where the w-space kernels Kpp, Kpo, Kgf and Kqq, given by the inverse Laplace transforms of the corresponding 
fc's, describe the smearing and growth of the original distributions Fs^v^Q"^) and G(u, Qg) through QCD radiation 
and splitting processes. 

The inverse Laplace transforms needed to implement Eqs. ([S]) and © can be calculated efficiently using the very 
accurate and extremely fast algorithms discussed in [1, these were used in the calculations reported here, and 
the results then converted to distributions in x and Q^. The numerical techniques needed are discussed in detail in 
the Appendix in [l[ . These allow the fast solution of the complete set of DGLAP evolution equations on a standard 
desktop or laptop computer. The kernel technique will be discussed elsewhere. 

The one-step inversion in Eqs. ([5]) and ^ is particularly useful in the case of devolution from large to small Q^: 
the variable r is then negative, the integrals that define Kpp and Kgg do not converge as ordinary integrals, and 
those kernels must be defined as generalized functions. This problem does not appear with the forms in Eqs. ([5]) and 
([6]) provided Fs{v, Q"^) and G{v, Q^) vanish sufficiently rapidly for — >■ that the products in Eqs. ^ and (jS]) vanish 
as a power of 1/s for s — > oo. These conditions are satisfied in practice. 

The evolved Fs(v, Q^) and G{v,Q^) must be continuous at quark thresholds where Uf changes. We treat the 
thresholds in as in |16l - [l8| . In the course of the evolution from the initial Qq to a larger final virtuality, may 
cross a threshold at = Mf where quark i becomes active, and the number nj oi active quarks increases by 1. 
This changes n/-dependent coefficients in the evolution equations. However, the continuity of Fs{v, Q"^) and G{v, Q^) 
as functions of is guaranteed if we evolve first from Qq to Mf, take the results at Q^ = Mf as new starting 
distributions, and then continue the evolution from Mf to Q^ with ri/ -> + 1. We otherwise neglect mass effects 
on the evolution. The same remarks apply to the case of devolution from Qf^ to a smaller Q^, with Uf then decreasing 
by 1 at each transition. 

We have checked that our methods accurately reproduce the LO results of CTEQ6L [l^ for the evolution of Fs and 
G when we use starting distributions taken from their published results. The errors in the evolved distributions are 
< 0.05% for CTEQ6L, as discussed in [J. Similarly, we reproduce the results of MSTW2008LO [13 for the evolved 
F, and G to < 0.1 - 0.5% Q. 

The solution of the nonsinglet evolution equations for quark distributions such as xq^ {x,Q^) = 
X [qi{x,Q^) — qi{x,Q^)] is simpler because of the absence of any coupling to the gluon distribution. The results 
in LO have the form 

Fns{v,Q^) = [fc„s(.s,r)/„,,o(s);i'] , (9) 

where fc„s(s,T) is the common LO non singlet evolution kernel and fns.o ~ ^ -fns(v,Qo)i 

We have discussed the generalization of these results to next-to- leading order in [2| . The decoupling of the evolution 
equations in that case requires a double Laplace transform and is considerably more complicated in detail, but can 
still be carried through analytically. We will not pursue that here. 
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B. Determination of the initial distributions 

In the following sections, we will apply our methods to an analysis of the combined HERA data [l^ on deep 
inelastic ep scattering. Those data determine the behavior of F^^(a;, Q^) very well for x < 0.1 for a wide range of 
Q^. F2^{x, Q"^) can therefore be taken as accurately known throughout that region through a global fit to the HERA 
data. 

Our objective is to check the consistency of LO QCD evolution with the HERA data by starting at an initial 
Qg, and evolving or devolving to the final values of where we can compare the evolved F2^{x,Q^) directly to 
the experimental results. To do this, we need to determine the initial gluon distribution G(x, Qq), which is not 
measured directly, and the initial singlet distribution Fs{x,Qq), both over the entire range {x, 1), evolve or devolve 
the distributions as discussed above, and then convert the resulting Fs{x,Q^) back to F2^{x,Q^). We will discuss 
the elements of this procedure in the following subsections. 



Determination of Go{x) — G{x,Qo) 



The LO evolution equation for F^^ix., Q^), easily constructed from the evolution equations for the individual quark 
distributions and the relation F2^{x, Q^) = J2i {n + 9*) {^^ Q^)j is 

Att dF^P{x,Q^) 



a,(Q2) dln{Q 



16 , 



(10) 



We have shown elsewhere [lj,ll5| that, assuming that a LO treatment of the DGLAP evolution of is consistent, 
we can invert Eq. pU)) to obtain G{x, Q^) at any given x, directly from a global fit to F2^{x, Q^) which includes 
the interval (x, 1) and a range of around the desired value. In particular. 



G(x,g2) =3 J'J'(x,Q2) 



dTT{x,Q^) 
dx 



- / FF 



3/2 ( 6 

7f 



sm 



%/7. z 
— — in — 

2 X 



2 cos 



V7, z 
— — in — 

2 X 



(11) 



where FF{x, Q^) is the function 



FF{x,Q'') 



4^ 



a,(Q2) ain(Q2) 



-AF^{x,Q' 



16 dF:'p 



dz 



(z,Q2)in 



dz + —X 
3 



(12) 



obtained by combining all the Fj'^-dependent terms in Eq. pO| and dividing the result by ^1- 

Since FF is determined by i^^J^, Eq. (|lip determines G directly from experiment provided the assumption of LO 
evolution is valid. We have found that the result for G{x,Q'^) at small x is fairly insensitive to the behavior of 
F2^{x,Q'^) at large x, so G{x,Q^) is determined at small x primarily by the HERA data. However, to get precise 
results, we need a global fit to F2^ that extends to x = 1. We will discuss that extension below. 

If LO evolution is consistent with the HERA data, the distribution G{x, Qq) determined by Eq. ([TT|) should satisfy 
the gluon evolution equation. We observed very early in our analysis that this condition was not satisfied. In 
particular, the derivative dG{x,Q^)/d\nQ^ was not equal to the sum of terms in the gluon evolution equation that 
involve weighted integrals of G and Fg. While this indicated that the assumption of LO evolution was not consistent, 
the strength of this conclusion was limited by the somewhat-limited accuracy with which the derivative of G could 
be determined. We have therefore adopted the alternative approach that we pursue here, and limit our consistency 
tests to the evolution of i^^J^, where direct comparisons with the HERA data are possible. 



2. Determination of the singlet distribution Fs{x, 



In the LO CTEQ6L [T6| and MSTW2008LO [T3l analyses which we will use for comparisons, the singlet quark 
distribution function Fs{x,Q^) was determined through a simultaneous fit to all the quark distributions and the 
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(18) 

), nf=4, (19) 



gluon distribution. Those analyses used earlier variations of the HERA data [lOKlJl '^^ combination with other 
data on deep inelastic electron and neutrino scattering that arc concentrated at higher x. Because of apparent 
incompatibilities among various data sets discussed in [3, [i3| , and the high accuracy of the combined HERA data at 
small X, we will adopt instead a hybrid approach in which we write Fs{x,Q^) in terms of F2^[x,Q'^) and relatively 
small nonsinglet quark distributions. We will then take F2^{x, Q^) from a global fit to the combined HERA data, and 
will use the nonsinglet contributions obtained in the CTEQ6L and MSTW2008LO analyses to construct Fs{x,Q'^). 
Those analyses differ in their treatments of as in NLO and LO, respectively. 
Introducing the nonsinglet quark distributions [l^ 

= x{q,-q,), i= 1,2,3, (13) 

T3 = x{u + u- d - d), (14) 

Tg = x{u + u + d + d-2s-2s), (15) 

lis = a;(u + It + d + (J + s + s — 3c — 3c), (16) 

T'24 = x{u + u + d + d + s + s + c + c- Ab- W), (17) 

we can write Fs{x, Q^) for different numbers n/ of active quarks as 

FAx,Q') = -^F^P{x,Q')-^Ux,Q^)~-Ux,Q^) + -T,,{x,Q' 
Fs{x,Q') = ^F^P{x,Q')-^Ux,Q')-^Ts{x,Q^) + ^T,,{x,Q')-^T2i{x,Q^), Uf^b. (20) 

Our procedure is now the following. We start with our global fit to F2^{x,Q'^) and pick an initial value of in 
a region where is well determined, here Qq = 4.5 GeV^, a value between the charm and bottom thresholds. We 
start by using the nonsinglet distributions T3, Tg, and T15 from the CTEQ6L (or MSTW2008LO) fit to the older 
HERA and high-x data to get an initial result for the singlet distribution Fs{x,Qq) from F2^{x,Q'^) using Eq. (|19p . 
We also determine G{x, Ql) from F2^{x, Qq) using Eq. ([TT|) . 

We next devolve Fs{x, Q^) to the charm quark threshold at = M^. The c and c distributions should vanish at 
(52 = M^, with Ti5(a;,M2) = Fs{x,Ml) for n/ = 3. However, because we have started with somewhat different data 
on F2^{x, Q^) than used in earlier analyses, this threshold condition will not be satisfied exactly. We therefore use the 
continuity of Fs{x^ Q^) at the nf = 3, Uf = A transition, set Ti5(a;, M^) equal to the devolved Fs{x, M^) for n/ = 4, 
and evolve T15 back to Qo using the LO nonsinglet procedure discussed in 2| to obtain a modified Ti^^x, Qg). This is 
used to get a modified Fs{x, Qq) from F2^{x, Qq), and the process is repeated if necessary until the result for Ti.g does 
not change significantly. The changes in Fg introduced by this procedure are small except near the charm threshold. 
The Fs{x,Ql) obtained from F2!^{x,Ql) using the modified Tig and the CTEQ6L (or MSTW2008LO) distributions 
T^^XjQq) and Ts{x,Qq) gives the initial singlet distribution for use in our subsequent calculations. 

The situation with respect to T24{x, Q^) is simpler. This distribution comes in at the bb threshold, where T24{M^) ~ 
Fs{x,M^) for Uf = 4. We therefore determine the initial distribution T24{x,M^) by evolving Fs(a;, Q^) from Qq to 
M^, and its extension to higher Q^, by evolving T24 from to using the results of [2| restricted to LO for 
nonsinglet evolution. 

Finally, the evolved or devolved F2^{x, Q^) is determined from evolved or devolved Fs{x, Q^) for a given n/ using 
the appropriate one of Eqs. (|18l) - (pn)) . The corresponding quark distributions can be obtained from Fs{x, Q^) and the 
(modified) nonsinglet distributions, as discussed later. 



C. A global fit to the combined HERA data for F^"'{x,Q'^) 



The constructions above depend on our having a global fit to the x and dependence of the structure function 
F^^{x, Q^). Berger, Block and Tan showed that ZEUS data from HERA [13, Ei could be parametrized accurately 
as a function of x and for a; < 0.1 by an expression of the form 



Fi{x,Q')^{l-x) 



Fp 



A\-a 



xp{l — x) 



Bin' 



Xp(l — x) 



1 — xp x{l — Xp) x{l — Xp) 

We will use the same parametrization for the complete HERA data sets as combined in [l 



(21) 
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In the expression in Eq. (|21|) , xp specifies the location in x of an approximate fixed point observed in the data where 
curves of F^p{x, Q^) for different cross. At that point, dF^P{xp, Q^)/d\nQ^ « for all Q^. = F^p{xp, Q^) is 
the common value of F2^. The dependence of Q^) is given in those fits by 

A{Q^) = ao + ai\nQ^ + a2\n^Q^, B{Q^) = bo + bi \nQ^ + b2 In^ ■ (22) 

We used this parametrization to fit the combined HERA data for > 1 GeV^. These data included 34 different 
Q2 values with x < 0.11, specifically, = 0.85, 1.2, 1.5, 2.0, 2.7, 3.5, 4.5, 6.5, 8.5, 10, 12, 15, 18, 22, 27, 35, 45, 60, 
70, 90, 120, 150, 200, 250, 300, 400, 500, 650, 800, 1000, 1200, 1500, 2000, and 3000 GeV^. The scaling point value 
Xp = 0.11 was taken to be fixed. 

The data set has a total of 356 datum points. The use of the sieve algorithm to sift the data to eliminate outliers 
as described in [2lj eliminated 14 points whose contribution to the of the fit was 125.0, roughly a quarter of the 
total. The values of the 7 fit parameters, along with their statistical errors, are given in Table HI The fit using the 
sieve algorithm gives a minimum with Xmin = 352.8. This must be corrected by the sieve factor TZ = 1.109 to account 
for the change in normalization of the function [2l|. This gives a corrected value TZ x Xmin = 391.4, so a corrected 
per degree of freedom of 1.17, a reasonable result for this much data. 

For the 296 points with > 2.7 GcV^ that we will consider later, the fit is excellent, with = 295. For 
comparison, the CTEQ6L ^ and MSTW2008LO [13 fits, made using the separate HI ^ and ZEUS [13, [uI data 
rather than the combined results, give of 3339 and 1329, respectively, with uncorrected values of x^/d.o.f. of 11.3 
and 4.5. 

Curves of the fitted F^^(a;, Q^) plotted as a function x are compared with the data in Fig. [T]for 24 values of Q^. 
The quality of the fit is evident. 

We emphasize that our fitting procedure is quite different from that used in other analyses. Our fit is directly to 
F2^{x,Q^), and its adequacy can be tested at that level. An investigation of possible alternative models with more 
parameters gave essentially equivalent results in the experimental region. We use the model in Eq. (|2ip and Eq. 
because of its simplicity and its reasonable behavior for small x [20| , and more importantly, for its excellent fit with 
a minimum number of parameters. With this approach, our fit to the HERA data is independent of any assumptions 
about QCD evolution, and will allow us later to obtain a direct test of the validity of purely LO evolution. In contrast, 
the usual methods, such as those in [iB - flst . start by assuming the validity of QCD evolution to some order in the 
strong coupling a.,, calculate -F^^ from a complete set of parton distributions evolved from some initial Qq, and then 
attempt to fit the data by adjusting the (many) parameters in the initial parton distributions. 

TABLE I: Results of a 7-parameter fit to the HERA combined data for F^^{x,Q^) for 0.85 < < 3000 GeV^. 



Parameters 


Values 




ao 


-8.471 X 10"^ ±2.62 x 10" 


-3 


ai 


4.190 X 10"^ ± 1.56 X 10' 


-3 


02 


-3.976 X 10"^ ±2.13 x 10" 


-4 


bo 


1.292 X 10"^ ±3.62 x 10" 


-4 


bi 


2.473 X 10"" ± 2.46 x 10" 


-4 


b2 


1.642 X 10"-' ±5.52 x 10" 


-5 


Fp 


0.413 ±0.003 


Xmin 


352.8 




7^ X Xmin 


391.4 




d.o.f. 


335 




7^ X Xmin/d.O.f. 


1.17 



D. Extension of the fit to high x 

Our fit to the data on F2^{x, Q'^) is so far restricted to the region x < xp; we have not attempted to fit the DIS 
data for x > xp from other experiments. Since the expressions for the evolved Fs{x,Q^) and G{x,Q^) in terms of 
their initial distributions at Qq given in Eqs. ([7]) and (|8]), and that for G in terms of given in Eq. ([TT|) . involve 
integrals that extend to a; = 1, we need also to extend the parametrization of F2^{x, Q^) to a; = 1. We will again use 
the results of earlier analyses, this time less directly, in making the extension. 



7 




FIG. 1: Comparison of our fit to the proton structure function F^^(a;,Q'^) with the combined HERA data [l^, plotted 
as functions of the Bjorken variable x, with increasing from the bottom to the top curves in each panel: (a) = 
3.5, 6.5, 15, 27, 120, 650 GeV^ (b) = 4.5, 10, 22, 45, 150, 800 GeV^; (c) = 0.85, 2.7, 12, 35, 90, 400 GeV^ (d) = 
1.5, 8.5, 18, 70, 250, 1200 GeV^. The fixed point in the fit was taken as xp = 0.11. 

We have found that the CTEQ6L and MSTW2008LO versions of F^^{x, Ql) for Ql = 4.5 GeV^ are well approxi- 
mated at large x by expressions of the form 

\xq J \1 - xq J 1 + axQ + oxq 

We will use this form to extend our fit to F2^{x, Q^) to the high-x region, where the HERA data are restricted to 
values of much larger than our chosen Qq, and F^^^x^Qq) is not well determined. In making this extension, we 
must choose the starting Xq sufficiently small that we avoid problems with our lack of precise knowledge of the x and 
dependence of F2^{x, Q^) for x near the fixed point in our fit. We have used xq = 0.03 in the present calculations. 
With this choice, the CTEQ6L result for F^'' is well fitted with a = 6.83, b = 13.0, and n = 3.75 in Eq. For 
MSTW2008LO, a = 4.83, b = 13.7, and n = 3.66. 

We match the expression in Eq. (|23|) in value and slope at xq = 0.03 to the expression in Eq. pT|) which describes the 
HERA data by adjusting the parameters Fq and ^, retaining the initial values of a, 5, and n. The changes necessary 
in jj, are fairly small, with increases of 4.6% and 4.0% in magnitude from the CTEQ6L and MSTW2008LO values, 
respectively. The changes in the normalizations are somewhat larger, 11.8% and 8.7%. To a good approximation, 
the extended distributions in the region x > 0.03 are simply scalings of the CTEQ6L and MSTW2008LO results for 
F2^{x,Qq), retaining the shapes of those distributions. Our final results at small x are insensitive to the details of 
these extensions. 

The determination of the initial gluon distribution at Qo involves further complications. As discussed in Sec. lIIBTl 
G{x,Qq) can be determined directly from F2'^{x,Q^). It can be shown from Eqs. pT|) and ([12]) that G{x,Qq) is 
actually determined mainly by F2^{x, Qq) and its derivative dF2^{x, Q^)/d\nQ^ at Qq] the integral terms in Eq. (fT2|) 
are small. The need to know dF2^{x, Q^)/d\nQ^ introduces some complication because the fixed point imposed in 
Eq. (|21|) reflects the observed dependence of F2^(x, Q^) for x near xp = 0.11 only qualitatively, and not precisely. 
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The HERA data near that point are restricted to >> Qq, and do not determine (5^)/91nQ^ in the region 

~ Qo where it is needed. The derivative at Qq is, in fact, only determined well by the fit to the HERA data for 
X << xp and ^ Qq- In particular, the fit to F^'' and its extension to high x do not give reliable results on its 
dependence for x >> 0.03. As a result, the expression in Eq. ([TT|) cannot be used to determine G in that region. 

We therefore adopt an approach similar to that used with F^^. We choose a small value of xq, xq = 0.03 where 
and its dependence are well determined, and determine G{x, Qo) for x < xq from the fit to F2^{x, Q^) using 
Eq. (jlip . The small uncertainties in the extensions of F^'' to large x affect only the integral terms in Eqs. (fTTj) and 
P^ . and do not affect the result for G significantly in the region of concern, x < 0.03. 

To extend the result for G to higher x, we fit the shapes of the gluon distributions G{x, Qq) given by CTEQ6L 
and MSTW2008LO for x > xq = 0.03 using the same functional form as in Eq. (|23|) . We use the results to extend 
G to a; > Xq by adjusting the analogs of the parameters fj, and Fq so that the extensions match the G derived for 
X < Xq in magnitude and slope at x = xq. The result is a gluon distribution G{x, Qq) that retains the basic shape of 
the CTEQ6L or MSTW2008LO gluon distribution for x > xq, merges smoothly into the form derived from F^^^ for 
X < Xq, and, in contrast to other analyses, involves no a priori assumptions about the form of G in the latter region. 



III. APPLICATIONS TO THE HERA DATA ON F^''{x,Q^) 

In this section, we summarize the results we obtained by applying our methods to an analysis of the HERA data 
on deep inelastic electron-proton scattering as combined by the HI and ZEUS experimental groups [T3| . 

We first examine the consistency of our results for Q^), G{x, Q^) , and the quark distributions with other LO 

results, represented here by CTEQ6L and MSTW2008LO. We find qualitative, but not quantitative agreement, with 
our evolved F2''{x,Q^) agreeing much better with the combined HERA data, and our G{x,Q-^) generally increasing 
much less rapidly at small x than the distributions found elsewhere. These changes will affect the results of cross 
section and other calculations performed using LO quark and gluon distributions. 

We then turn to a central question, the consistency of a LO treatment of the QCD evolution, and examine the 
consistency of the structure function F^^(a;, Q^) determined by LO evolution with the HERA data. We conclude on 
the basis of systematic, Q^ -dependent discrepancies that LO evolution cannot give an adequate description of those 
data. At least NLO corrections are needed. We emphasize that this conclusion is independent of any calculation of 
the NLO corrections, and follows directly from the Q^ dependence of the data. 



A. Basic results and comparisons with other analyses 

1. Starting distributions and sum-rule tests 

Our results are based on the smooth global fit to the measured Q^) discussed in Sec. Ill CI The fit was very 

good, as seen in Fig. [1] and determined our starting distributions at Qq ~ 4.5 GeV^, a value chosen in the region of 
dense data where the x and Q^ dependence of are tightly constrained. 

F^^{x, Ql) is fixed by the fit. We determin ed th e initial G(x, Qg) directly from the fit to F^^{x, Q"^) using Eq. (fTTj) 
and the extensions to high x discussed in Sec. | II Dl The uncertainty in our derived G at small x is determined mainly 
by aF2^^'(x,Q2)/aing2, and is quite small [13" We compare these initial distributions with those that resulted from 
the CTEQ6L and MSTW2008LO analyses in Fig. [2l There are clearly significant differences in the magnitudes and 
X dependence of distributions among the sets. 

We note first in Fig. [2{a) that the difference between the extensions of F^^(a;, Qq) for x > 0.03 we obtain for 
CTEQ6L-like and MSTW2008LO-like shapes is very small. These differences lead to negligible effects in the evolution 
of and G at small x. The differences evident between our curve for and those shown for CTEQ and MSTW 
in Fig. [D^a) result from their failure to fit this quantity accurately, presumably attributable in part to their use of the 
older HI and ZEUS versions of the data. 

The differences in our curves for G in Fig.[2l[b) from those of the CTEQ6L and MSTW2008LO analyses result from 
the difference between their F^^ and ours. The marked difference between the curves shown for our CTEQ-like and 
MSTW-like gluon distributions results from the different treatments of as used by the two groups, which we follow 
here. CTEQ6L treats to NLO, with 



/3oln(QVA2) 



2/3i ln[ln(QVA^)] 
Pi ln(QVA2) 



(24) 



2 19 
/3o(n/) = 11-3"/, /3i(n/) = 5I-yn/. (25) 
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The value of as is fixed to the measured value at the Z-boson mass, as{M^) = 0.118 for n/ = 5, and the value of 
A(n/) is then adjusted at the b and c thresholds where Uf decreases by 1 to assure continuity. 

MSTW2008LO, in contrast, uses only the first, LO, term in Eq. (|25|) for presumed consistency in a LO analysis, and 
treats the value of as at = 1 GeV^ as a parameter in their fitting procedure. This leads to a value as(Mf ) = 0.139. 
The two versions of as do not agree well, with the MSTW2008LO version b eing considerable larger at all Q^. We 
note that the dependence of as is actually well determined by experiment [22|, with the results well described by 
the NLO expression [2^ for as fixed to as{M^). Since is also known, the CTEQ-like determination of G{x,Qq) 
is based entirely on measured quantities, with the assumption that a LO analysis of the evolution is adequate. Our 
MSTW-like approach uses the MSTW2008LO version of as , but at the expense of poor agreement with the measured 
dependence of a^. 

Figures l^c) and (d) show the extensions of the curves in 2(a) and 2(b) to small x. We emphasize that with the 
assumption that the LO evolution equation for i^^'' is satisfied, a necessary condition for a consistent LO analysis, 
our initial gluon distribution G at Qg = 4.5 GeV^ follows directly from our global fit to the x and dependence 
of the HERA data on F2^{x,Q^) and its extension to large x. In this sense, F2^{x,Ql), G{x,Qq), and up to small 
corrections, Fs{x,Qq) are all determined by experiment for x > 10~^ where there are substantial HERA data, and 
determined to lesser accuracy down to a; ~ 10~^ where the data at presumably perturbative values of run out. It 
is not necessary to determine these quantities indirectly through initial parametrizations of the complete set of quark 
distributions and G, with the many parameters determined only in a fit to the data. 

We conclude that the strong divergences of and G evident in the MSTW2008LO curves in Figs, djc) and (d) 
are not realistic in a LO analysis. The lesser differences between the CTEQ6L results and ours in Figs. [IJc) and (d) 
are mainly in the region x < 10~^ where some extrapolation from the data is necessary, so it is less definitive. 




FIG. 2: Comparison of our starting distributions F^^{x,Qo) and G{x,Qo) with those of CTEQ6L and MSTW2008LO at 
Qo = 4.5 GeV^. (a) F2^ from our fit to the HERA data, extended to a; > 0.03 using the method described in the text based on 
the CTEQ-like (red dashed lines) or MSTW-like (blue dotted lines) shape of F^^ at larger x. The original CTEQ6L (red solid 
lines) and MSTW2008LO (blue dot-dashed lines) versions of F^'' are shown for comparison, (b) The G{x,Ql) derived from 
F2^ using the condition that F^'' satisfy its DGLAP evolution equation in LO, using the NLO (red dashed lines) or LO (blue 
short dashed) versions of Qs, compared to the corresponding CTEQ6L (red solid lines) and MSTW2008LO (blue dot-dashed 
lines) distributions, (c) Extension of (a) to small x. The CTEQ-like and MSTW-like shapes for fj''' for x > 0.03 cannot be 
distinguished on the scale of the figure, and only the former is shown, (d) Extension of (b) to small x. 
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Following the procedures discussed in Sec. Ill B"2l we used the fit to Q^) and the results for the nonsinglet 

quark distributions Vi, T3, Tg and the initial given by CTEQ6L or MSTW2008LO, to determine the corresponding 
LO result for Fs{x,Ql). 

As a test of our procedures, we evaluated the QCD momentum sum rule, which should give 

/ dx[F4x,Ql) + G{x,Ql)] =1. (26) 
Jo 

We find that it is satisfied to ^ 0.1% (1.2%) at Q"^ = 4.5 GeV^ for the Fg and G derived from the extended fit to 
F:P' using the nonsinglet distributions from CTEQ6L (MSTW2008LO) and the method of Sec. Ill B 21 Because of 
the structure of the splitting functions, the sum rule for the evolved distributions is automatically satisfied to similar 
accuracy at all . 

This resuh may seem startling: the CTEQ6L and MSTW2008LO resuhs for i^s(a;, Qg) and G{x, Ql) also satisfy the 
sum rule, used in those fits as a constraint, but our global fit to the combined HERA data at Qq lies considerably above 
the F;^^(a;, Qq) calculated from their quark distributions as seen in Fig. [2l with similar differences in Fg. However, 
the gluon distribution G{x, Qq), calculated from the requirement that F^^ satisfy its LO DGLAP evolution equation 
exactly, is smaller than the G obtained in other analyses in the region of x that contributes significantly to the sum 
rule, as seen in Fig. [21 

The two effects compensate for each other numerically. The contributions to the momentum sum rule from Fs and 
G &t Ql = 4.5 GeV^ are 0.550 (0.622) and 0.452 (0.377) for the CTEQ6L (CTEQ6L-like) distributions, with the 
calculated sum rule equal to 1.002 (0.999). The resuhs for the MSTW2008LO (MSTW-like) distributions are similar, 
with contributions to the sum rule from F^ and G of 0.565 (0.637) and 0.434 (0.373) at Ql = 4.5 GeV^, for total of 
0.999 (1.010). We have not used the sum rule as a constraint, as is done in other analyses. Its satisfaction follows 
from the data and our determination of G in terms of F2^ . We conclude that our extensions of F^^ and G to the 
large-x region cause no problems. 

The quark number sum rules 

dx{u-u){x,Ql) ^2, I dx{d-d){x,Ql)^l, (27) 

Jq 

are different. Because we set the nonsinglet distributions u — u and d — d equal to the corresponding CTEQ6L or 
MSTW2008LO distributions and do not change them in our hybrid analysis, the quark number sum rules are satisfied 
automatically to the extent that they were satisfied by the CTEQ6L and MSTW distributions, namely to ^ 0.4% 
(~ 0.5%). The changes introduced in the separate u and u, and d and d distributions by the changes in F2'^ and Fg, 
are confined to the singlet combinations u + u and d + d, and cancel in the differences u — u and d — d. 

The corresponding sum rules for s — s, c — c, and b—h give zero in the CTEQ6L-based analysis since those quarks are 
produced only in pairs through gluon splitting. For the MSTW2008LO-based input, s ^ s initially. The very small 
difference is not changed in our analysis because we keep the nonsinglet distributions fixed, and the strange-quark 
sum rule remains constant at w 0.0028. The c, c and b, b quarks are produced only in pairs, and the quark sum rules 
give zero. 



2. Leading-order gluon and quark distributions 

We evolved the starting distributions for Fg and G from Ql to lower and higher values of Q^ using the Laplace 
transform methods sketched in Sec. Ill Al using the numerical techniques discussed in the Appendix to [§]. We compare 
the evolved gluon distributions G{x,Q^) to those of CTEQ6L and MSTW2008LO ^ in Fig.El 



It is evident from the figure that our gluon distributions are somewhat smaller than those of CTEQ6L and 
MSTW2008LO, quite significantly so for the latter at small values of x where MSTW uses a strongly power-law 
divergent parametrization with their initial G{x, Ql). Our CTEQ- and MSTW- based results also differ significantly, 
the result of the differing initial distributions seen in Fig. [2] and the different treatments of as NLO and LO 
respectively. 

It is straightforward to combine our results for Fs{x,Q^) with the original nonsinglet distributions Vi, T^, Tg, and 
the modified Tis and T24, Eqs. p^ -lfTT ]) to obtain the quark distributions that lead to these results. The results differ 
from the individual quark distributions given by CTEQ6L and MSTW2008LO because of the changes in the HERA 
data, and, more importantly, because of our treatment of the starting distributions for the evolution of Fs{x, Q^) and 
G(x,Q2). 

Our results for the quark distributions are shown in Figs. U and [5] for the treatments based on the CTEQ6L and 
MSTW2008LO nonsinglet terms, respectively. The differences from the input distributions arc not large in the region 
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FIG. 3: Comparison of our evolved gluon distributions G{x,Q'^) with the CTEQ6L and MSTW2008LO distributions, (a) Our 
G (dashed curves) at — 10 (black), 35 (red), and 120 (blue) GeV'^ , bottom to top, compared to the CTEQ6L G (solid 
curves), (b) Our G (dotted curves) at = 10 (black), 35 (red), and 120 (blue) GeV^ , bottom to top, compared to the 
MSTW2008LO G (solid curves). 

of the HERA data, but some changes are evident at higher values of x, and, especially for MSTW, at very small x. We 
attribute the differences to the parametrizations of the quark and gluon distributions used by those authors, which 
have a strong power-law dependence on 1/x at small x, with the many parameters adjusted to fit the data used. 

Our method is based instead on our overall fit to F2^{x, Q^), and the information that can be derived from it. It 
uses the earlier nonsinglet distributions only in calculating small terms involved in the transitions between F^'' and 
Fg. The results on the fit shown in Fig. [1] suggest that its x and dependence are well determined for of a few 
GeV^ for x > 10~^. This allows the reliable derivation of the starting distributions needed in the solution of the LO 
evolution equations in that region. In that sense, our results are as reliable as allowed by the assumption of strict LO 
evolution. They do not depend on choices of parametrizations for initial quark and gluon distributions. The results 
shown in Figs. [31 HI and [5] follow. 

B. Check of the consistency of LO DGLAP evolution with the HERA data 

As a final application of our methods, we turn to the question of the consistency of LO evolution with experiment. 
We show that the structure functions F^^(a;, Q^) obtained by LO evolution from the initial distributions at Qq 
determined by the HERA data are not consistent with the data at higher and lower values oi Q^. A consistent 
analysis must therefore include higher-order terms in in the evolution equations, and distributions evolved out of 
the experimental region using the LO DGLAP equations cannot be used with confidence. 

We plot the ratios (^F^ evolved ^ ^2^hera) /^2^hera ^'^^ ^0 values of where there are data in Figs. [HI and [T] Here 
^2'ovoivod distribution evolved (or devolved) from Qq = 4.5 GeV^, and ^"^hera HERA data. 

We also show the ratios with -F2''hera replaced in the numerators by the actual data points. 

We can see from the figures that the evolved distributions differ systematically from the fit and the data, falling 
too low for > Qq for x in the range ^ 5 x 10~^ — 10~^, and too high for x < 5 x 10~^. The discrepancies increase 
systematically with increasing Q^, span about a 10% range for 0.001 < a; < 0.01, and have the same pattern for the 
analyses based on the CTEQ6L and MSTW2008LO nonsinglet distributions. The datum points follow the curves, as 
they should; the problem is not in the fit. The systematic increase of the discrepancies with increasing Q^ indicates 
that they are the result of incorrect evolution at LO, with the evolved F2^ not growing sufficiently rapidly with Q^. 
We conclude that LO DGLAP evolution of F2^ is inconsistent with the combined HERA data. 

The systematic trends are evident quantitatively in Table [III Using the 296 data points in our sample of the 
combined HERA data for Q^ > 2.7 GeV^, we find a (x^ P^r degree of freedom) of 295 (0.996) for our fit from Sec. 
Ill C[ 1480 (5.00) for the evolved that used the CTEQ6L nonsinglet terms to convert between and Fg, and 
502 (1.70) for the evolved F^^ that used the nonsinglet distributions of MSTW2008LO. Our direct fit to the HERA 
data is quite good given the large amount of data, with probability P = 0.126 when is corrected for the sieve factor 
[2H TZ = 1.109. The evolved distributions have essentially zero probabilities of being correct statistically. 

The difference in the values of for the CTEQ6L- and MSTW2008LO-based treatments of the nonsinglet terms 
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FIG. 4: Plots of the quark distributions obtained by our method using the nonsinglet distributions from CTEQ6L [T^, shown 
for: (a), = = 1.69 GeV^ (b), = 10 GeV^ (c), = 35 GeV^ and (d), = 120 GeV^ The sohd hnes give our 
distributions. The dashed hnes are the CTEQ6L distributions. The curves show xq{x,Q^) for, top to bottom in each panel at 
X ~ 0.1, the u (red), u (blue), d (black), d (green), s = s (purple), c = c (orange), and b = b (brown) quarks 



is the result primarily of the different treatments of in the two cases. The NLO treatment in CTEQ6L is fixed to 
the value of as at M|, and agrees well with the measured values of ag down to M^- In contrast, the value of the LO 
version of as at = 1 GeV^ is used in MSTW2008LO as a fitting parameter. The result is an as that is larger than 
the NLO version by about 40% at = 1 GeV^, and 18% at M|, so it does not agree with the measured values. This 
results in rather different starting distributions at Qq = 4.5 GeV^ in the two cases, as seen in Fig. [2l and to more 
rapid QCD evolution in the case of the MSTW2008LO-based treatment. Although the resulting is reduced, the 
systematic problems with the evolved remain, as seen in Fig. [71 and the result is still unacceptable statistically. 

This failure of LO evolution to give an accurate description of the separate HI and ZEUS data has been noted in 
[13, [13 J s-nd no doubt elsewhere, in connection with poor values of the for obtained in LO in those analyses, 
and the improvements afforded by a NLO treatment of the parton distributions. The systematic nature of the problem 
is somewhat obscured there by the way initial conditions are imposed through many-parameter descriptions of the 
complete set of parton distributions, and the subsequent adjustment of those parameters to minimize the of the 
fit. 



IV. SUMMARY AND CONCLUSIONS 



In the present paper, we have applied recently developed methods based on Laplace transforms to a LO analysis of 
the HERA data on deep inelastic ep scattering as combined by the HI and ZEUS experimental groups [l^. We have 
used a hybrid method, in which we convert the measured structure function F2 ^{x,Q'^) to the singlet distribution 
Fs{x, Q^) which enters the evolution equations, taking the small contributions of nonsinglet quark distributions to this 
conversion from other analyses, and extending the fit to the HERA data for x < 0.1 to .t = 1 using the shape of F2^ 
determined in those analyses. Here we used the results of the CTEQ6L and MSTW2008LO 117^ analyses, which 
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FIG. 5: Plots of the quark distributions obtained by our method using the nonsinglet distributions from MSTW2008LO [13], 
shown for: (a), = = 1.96 GeV^ (b), = 10 GeV^ (c), = 35 GeV^; and (d), = 120 GeV^ The sohd hnes give 
our distributions. The dashed lines are the MSTW2008LO distributions. The curves show xq{x, Q^) for, top to bottom in each 
panel at a; ~ 0.1, the u (red), u (blue), d (black), d (green), a (purple), a (magenta), c = c (orange), and 6 = 6 (brown) quarks. 



used the older HI [T^ and ZEUS [T3, [Hi data along with data from other experiments, mostly at higher values of x 
than the HERA data. This procedure determines the starting distribution Fs{x,Q'^) at the starting point Qq = 4.5 
GeV^ chosen for the DGLAP evolution. 

As shown earlier [13, [iBl, the necessary starting distribution Gq{x) = Gq{x,Q})) for the coupled evolution of Fg 
and G can be obtained in LO directly from a global fit to the structure function F2^{x,Q'^) by requiring that the 
LO evolution equation for F2^{x,Q^) be satisfied for = Qq. Both F2o{x) = F2^{x,Qq) and Go{x) are therefore 
determined directly by experiment through our fit to the HERA data for x < 0.1 and its extension to higher x, 
without the need for a solution of the complete set of coupled parton evolution equations or any assumptions about 
the functional form of G. Our results at small x are insensitive to the details of the extensions. 

We picked a starting value Qq = 4.5 GeV^ for the evolution which is well within the region of dense data. We then 
solved the LO evolution equations using very fast and accurate methods discussed elsewhere d, and combined the 
evolved F^ with the evolved nonsinglet distributions of CTEQ6L and MSTW2008LO to obtain a new set of quark 
distributions. These differ from the quark distributions obtained in those analyses because of our use of the combined 
HERA data rather than the original HI and ZEUS results, and our different determination of the starting distributions 
in Fs and G for the evolution. The differences in the quark distributions are significant in some regions. Our gluon 
distributions differ markedly from those of MSTW2008LO at small x as seen in Figs. [2]and|3l 

Finally, we compared the evolved structure function F2^{x, Q^) to the HERA data as a test of the consistency of 
LO DGLAP evolution. The initial distributions of Fs and G at Qq were determined by up to the small nonsinglet 
corrections, and were consistent with LO evolution by construction. We concluded that LO evolution is actually not 
consistent with those data on the basis of systematic trends evident in the evolved distributions. This conclusion 
docs not depend on the explicit calculation of NLO effects. It is supported by a analysis, but in contrast to other 
approaches, we could not attempt to reduce the by adjusting the shapes of the initial distributions: we had no 
arbitrary parameters to adjust. 
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TABLE II: The of the F^^' fr om our fit to the combined HERA data, and of the evolved F^'' obtained by LO evolution 
from Qq = 4.5 GeV^. The starting distribution G(x,Qo) was derived from the fit to F^^. The initial F^i.Qo) 

was obtained 

from F2^ using nonsinglet corrections from the CTEQ6L [l^ and MSTW2008LO [l3l analyses. The last lines give sums of all 
the rows above them except for the starting value Qq = 4.5 GeV^. 



Q (m (_.eV ) 


No. of datum points 


X (our fit) 


Xovolvcd) CTEQ 

corrections 


xLoivcd, MSTW 
corrections 


2.7 


9 


10.0 


15.4 


28.0 


3.5 


9 


11.1 


11.6 


11.5 


4.5 


9 


6.1 


6.1 


6.1 


6.5 


13 


14.2 


13.6 


14.3 


8.5 


9 


7.6 


7.6 


10.7 


10 


7 


2.4 


3.8 


7.1 


12 


10 


11.5 


15.1 


19.5 


15 


10 


10.6 


5.0 


28.5 


18 


9 


2.74 


25.0 


15.9 


22 


9 


12.4 


14.0 


10.4 


27 


12 


9.1 


52.9 


15.4 


35 


11 


8.8 


81.2 


11.0 


45 


11 


8.0 


96.9 


7.3 


60 


10 


17.2 


158.9 


19.4 


70 


9 


13.7 


68.7 


12.6 


90 


11 


13.0 


175.9 


49.1 


120 


12 


6.8 


102.0 


25.1 


150 


12 


15.9 


71.0 


16.5 


200 


14 


21.0 


114.6 


33.5 


250 


14 


15.8 


86.6 


25.6 


300 


15 


18.9 


83.5 


24.4 


400 


14 


18.7 


76.0 


21.6 


500 


11 


5.8 


37.7 


18.7 


650 


12 


10.4 


57.5 


21.4 


800 


9 


11.0 


40.8 


18.8 


1000 


9 


6.1 


13.3 


5.9 


1200 


9 


10.0 


33.1 


18.2 


1500 


6 


5.8 


9.9 


5.7 


2000 


5 


0.33 


0.26 


1.1 


3000 


5 


6.3 


7.5 


4.7 


Sum (without = 4.5) 
xVd.o.f. 


296 


295.2 
1.003 


1480 
5.00 


502 
1.70 



In the Appendix, we give an equally accurate, though approximate, method which works directly with the exact 
DGLAP LO evolution equation for i^^^ coupled to an approximate evolution equation for G. This approach is 
independent of the nonsinglet distributions, and its implementation uses only the experimental results as extended 
above. The results of the analysis are the same: LO evolution of F^^ is inconsistent with the HERA data. 
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Appendix A: Approximate coupled evolution equations for F2^{x,Q^) and G(x,Q^) 

In this Appendix, we point out that we can obtain a direct test of the adequacy of LO evolution using evolution 
equations coupling F2^\x,Q'^) and G{x,Q^). In particular, we use the exact LO evolution equation for F^^, and an 
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FIG. 6: Fractional accuracy ^J2'cvoivcd ~ /^2'hera evolved from Qq — 4.5 GeV^ relative to our fit to the combined 

HERA data [ij, compared to the same ratio with the data for F^'' used in the numerator. The initial and final conversions 
between and the singlet distribution Fs are based on the CTEQ6L nonsinglet quark distributions [l^. Results are given 
for (a) — 2.7 (black dots), 12 (red squares), 35 (green diamonds), 90 (blue triangles), 400 (orange inverted triangles) GeV^; 
(b) = 3.5 (black dots), 8.5 (red squares), 18 (green diamonds), 70 (blue triangles), 250 (orange inverted triangles) GeV^; (c) 
= 6.5 (black dots), 15 (red squares), 27 (green diamonds), 120 (blue triangles), 650 (orange inverted triangles) GeV^; and 
(d) = 10 (black dots), 22 (red squares), 45 ( green diamonds), 150 (blue triangles), 1200 (orange inverted triangles) GeV . 



approximate version of the evolution equation for G in which Fs is replaced by a multiple of F^''. 

The advantage of this approach is that it deals directly vifith the experimentally accessible function F^'^'lx, Q^), and 
gives a direct test of the adequacy of LO evolution with no input beyond a global fit to F^''. It does not require direct 
knowledge of the nonsinglet quark distributions, but correspondingly does not provide individual quark distributions 
unless Vi, Ts, Tg, T15, and T24 are known. If these are to be used, the method developed in Sec.|IT]is to be preferred. 

The results on the evolution of F2^{x, Q^) from its initial distribution at Qg = 4.5 GeV^ obtained by this method 
differ insignificantly from those obtained with the method in the body of the paper, with fractional differences small 
on the scale of the differences of the evolved F2^ from the data shown in Figs. |n]and[71 We conclude again that the 
assumption LO evolution is not consistent with the HERA data. 

We obtain our evolution equations for i^j'^^ and G as follows. The exact LO evolution equation for F^^(a;, Q^) is 
given in Eq. pUj) . This equation couples F2^ to the gluon distribution G. The exact evolution equation for G couples G 
instead to the singlet quark distribution Fs{x, Q"^) = x {qi + qi) {x, Q^), and not to F^'' = J2i ^i^ ili + Q^)- 
Fs{x,Q^) is not determined directly by experiment. However, we note that the nonsinglet contributions in the 
transition from to Fg given in Eqs. P^ - ipH)) are very small, and will simply replace Fs{x,Q'^) in the usual 
evolution equation for G{x, Q^) by the leading, F^^'-dependent terms in Eqs. p^ - ipU)) . Fs{x, Q"^) « a{nf))F^^{x, Q"^) 
with a(n/) = 18/5 for < Q"^ < M^, and 45/11 for < Q'^ < These relations are actually only expected 

to hold for well above thresholds, where the new quarks can be taken as fully excited; we will use them as stated. 

We use the resulting approximate evolution equation for G with the exact LO evolution equation for in Eq. (|10p , 
and solve for Fj''^ and G using the methods developed earlier [l|, H, [l3] . The accuracy of the method is evident from 
Table IIIIl where we compare the results for the evolved F^^ obtained using the approximate method with those 
obtained using the exact evolution equations for F, and G and the CTEQ6L nonsinglet corrections in the F^^^, Fg 
transition as described in Sec. Ill B 21 The accuracy is similar for the MSTW2008LO-based nonsinglet corrections. 
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FIG. 7: Fractional accuracy ~ -F^^hera ) /^2'heha of the structure function F,]^ {x,Q^) evolved from Qq = 4.5 GeV^ 

relative to our fit to the combined HERA data compared to the same ratio with the data for used in the numerator. 
The initial and final conversions between F^'' and the singlet distribution Fs are based on the MSTW2008LO nonsinglet quark 
distributions Results are given for (a) = 2.7 (black dots), 12 (red squares), 35 (green diamonds), 90 (blue triangles), 
400 (orange inverted triangles) GeV'^; (b) = 3.5 (black dots), 8.5 (red squares), 18 (green diamonds), 70 (blue triangles), 250 
(orange inverted triangles) GeV^; (c) = 6.5 (black dots), 15 (red squares), 27 (green diamonds), 120 (blue triangles), 650 
(orange inverted triangles) GeV'^; and (d) = 10 (black dots), 22 (red squares), 45 (green diamonds), 150 (blue triangles), 
1200 (orange inverted triangles) GeV^. 

The approximation of replacing Fs by a multiple of Fj^^ is only good to about 5-7% at = 5 GeV^, a value above 
the c-quark threshold but below the 5-quark threshold, and also at 100 GeV^, well above the b threshold, so the effect 
of these errors on the final is clearly greatly reduced by the nature of the evolution. We can understand this 
qualitatively as follows: the evolution of G at small x is driven mainly by G itself, which is accurately known at the 
initial Qq from the condition that the measured F2^{x,Qq) satisfy its evolution equation. The final errors in G are 
therefore small at small x, and their effect on F^^ is further suppressed by the contributions from F^^ itself to its 
evolution. In addition, G is small at large x, and errors in the approximate G in that region have little effect on the 
final F^^. Overall, the limited accuracy of the approximate Fg has only a small effect on the evolved G, and as a 
result, even less effect on the exact evolution of i^^^ from its known initial distribution. 

We have described the methods we use to solve the coupled evolution equations for Fs{v, Q^) and G(u, Q^) in detail 
elsewhere [H, [l]- We use the same methods here to solve the coupled equations for Q^) and G{v, Q^), so we only 
point out the changes. We begin with Eqs. ([2]) and ([3]) which express the Laplace transforms /2(s,t) and g{s,T) of 
the distribution functions in terms of their initial distributions and, here, a set of new kernels kij{s,T) — >■ k-^ {s,t). 

The kernels have the same form as those given in 0,0, but with the coefficient functions $, 6 that appear there 
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TABLE III: Fractional differences Af and Ac (in %) between the and G distributions obtained using the "exact" trans- 
formation between Fs and -F^^'' described in Sec. Ill B 21 and those obtained using the approximate method based on the exact 
evolution equation for -F^'', and an approximate gluon evolution equation with Fa replaced by a multiple of F^'' as described 
in this Appendix. The same starting distributions for F^'' and G at Qy = 4.5 GeV^ were used in both cases. In the "exact" 
method, we used the nonsinglet terms from CTEQ6L to convert between Fs and F^''. The last column shows the percentage 
rms differences between the distributions from the two methods for 10"^ < a: < 0.5. Results obtained using the nonsinglet 
terms from MSTW2008LO are very similar. 
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replaced by functions and 



= ^"3 (,7TT + ^ + ^^'^'^' + ^^ + ^^^j' ^^^^ 

Qf\s) ^ Yeff— + ^ V (A2) 

^ ^ ' ^ ' \s + l s + 2 s + 3 J ' ^ ^ 



1 1 



ef(.) - (---^ + -^). (A4) 



'2 2 1 

3"''"-" Vs ~ s + 1 ^7+2 

These functions differ from the corresponding functions in the case of F^, G in the coefficients in the 6's, hence the 

(2) 

introduction of the superscripts 2 to distinguish the two cases. The kernels k^j have the same formal structure as the 
original fc.y, and the final solutions are obtained as described in Sec. Ill A] using the very fast and accurate algorithms 
for calculating inverse Laplace transforms introduced in H, [^. The methods needed in practice are discussed in 
the Appendix of [l[. The results are essentially the same as those presented in Sec. IIIIBl and we draw the same 
conclusions as there. 
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